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We consider the radiation transfer problem in the discrete-ordinate, plane-parallel approach. We 
introduce two benchmark problems with exact known solutions and show that for strongly non- 
homogeneous media the homogeneous layers approximation can lead to errors of 10% in the estima- 
tion of the intensity. We propose and validate a general purpose numerical method that transforming 
the two-boundary problem into an initial boundary problem, using an adaptative step integration 
and an interpolation of the local optical properties, can improve the accuracy of the solution up to 
two orders of magnitude. This is furthermore of interest for practical applications, such as atmo- 
spheric radiation transfer, where the scattering and absorbing properties of the media vary strongly 
with height and are only known, based on measurements or models, at certain discrete points. 



I. INTRODUCTION 

The mathematical modeling of radiative transfer in which the phenomena of absorption, emission and scat- 
tering are taken into account is usuaUy made using the linearized Boltzmann equation, also known as radiation 
transfer equation. This equation describes the transfer of radiation, with a given wavelength, through a medium 
with certain absorbing and scattering properties. A particular application of this equation is the study of ra- 
diation transfer in the atmosphere, where the medium properties vary strongly with height. Getting accurate 
solutions of this problem is important for evaluating energy balance on planetary atmospheres as well as for 
the so called "inverse problem" where the boundary conditions (in particular the ground albedo properties) are 
deduced from measurements of radiation and knowledge of the medium characteristics. 

In this work we consider the time-independent, monochromatic radiative transfer equation using the well- 
tested and widely used discrete-ordinate method of Stamnes et al. HI and the plane parallel approach where 
the optical properties depend only on the vertical coordinate z. The procedure requires the solution of a system 
of n coupled linear ordinary differential equations (one for each stream or discrete-ordinate component of the 
intensity). This set of equations is subject to a two-point boundary condition at the top and bottom of the 
medium. In the general case no analytic solutions exist for this problem since the medium optical properties 
(phase function, absorption and scattering coefhcients) depend on the position z in the vertically inhomogeneous 
medium. 

To obtain a formal solution, the medium is generally assumed to be layered with piecewise constant optical 
properties, i.e. it is divided into N adjacent homogeneous layers where the eigenvalues are computed. The 
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FIG. 1: Layer with altitude dependant extinction and scattering coefficients. 

coefRcients of the solution are determined by imposing the continuity condition at the boundaries between 
adjacent layers and the two-point boundary conditions at the top and bottom of the medium |^|^. 

Unfortunately in many cases the optical properties of the medium are not homogeneous, in fact they show 
strong variations along the vertical axis and have therefore different characteristic length scales. In these cases, 
as we will see below, the homogeneous layers assumption may lead to errors of 10% in the estimation of the 
scattered radiation. Furthermore in most practical applications the medium characteristics are only known, 
based on measurements or calculations, at a discrete number of points 0|. It is therefore required to implement 
a method that is able to cope both with the strong variations on the medium properties and with the discrete 
character of the information available. 

Here we solve the discrete-ordinate approximation to the radiative transfer equation with a different approach. 
Since the equations and the boundary conditions are linear in the intensity, the two-point boundary problem 
can be solved with a shooting method. The problem is then reduced to the solution of an initial value problem 
which can be solved numerically. We propose a numerical integration of the initial value problem that applies 
an adaptive step method (such as step doubling), interpolates the optical properties from the discrete set of 
available data (for example using a cubic spline) and uses a weighted evaluation of the integrand along the 
interval (in our case a 5th order Runge-Kutta scheme). This allows the adaptation of the numerical integration 
to the local characteristic length-scales of the problem under consideration. The assumption of homogeneity is 
thus not required. 

As an example we apply this method to solve the two-stream discrete ordinate version of the, non-emitting, 
radiative transfer equation in one spatial dimension (these results can be equivalently extended to multi-stream 
and emitting versions of this equation) . We present two benchmark cases with known analytical solutions and 
compare the performance obtained when the step doubling interpolating method is used and when the piecewise 
homogeneity is imposed. We will show that using the same available discrete information, our method can 
significantly improve the accuracy of the solution. In section ^] we will first introduce the radiative transfer 
equation and the equations for the direct and diffuse intensity components in the two-stream discrete ordinate 
approach. In section ITTll we will describe two benchmark cases characterized by linear and exponential height 
dependence of the optical coefficients in the equations. In section Hvl we will explain our numerical procedure 
and will validate our results and those obtained under the assumption of piecewise homogeneous layers against 
the exact analytical solutions of the benchmark problems. Finally conclusions are presented in section IVI 



II. MULTIPLE SCATTERING AND THE RADIATIVE TRANSFER EQUATION 

Our aim is to solve the radiative transfer problem in the plane parallel approach. As represented in Figure ^ 
this corresponds to the case in which optical properties only depend on altitude z in a plane parallel geometry. 
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Following for instance the radiative transfer equation in this approach is written as, 

cose — = -(3T(z,X)I{z,e,(t),X) + — J d(j) (1) 



sine' de'p{z, 0, (j), z', 9', (/)', X)I{z, 0', (/)') 

where 9 and are respectively the polar and azimuthal angles. Here Psca represents the attenuation due to 
scattering effects and firiz, X) is the extinction coefficient defined as Pt{z, X) = Psca{z, X) + I3i{z, X) where 
the summation in i extents to all the molecular components considered. Here (3i{z,X) = ni{z)ai{X) is the 
absorption coefficient of a given component, ai{X) is the attenuation cross section due to absorption and ni{z) 
the atmospheric number density for species i (which generally depends exponentially on the height z). In 
Eq. f^l, Piz,d,4'j z' ,9' , (j)' , X) is the phase function of the scattering particles which is normalized as follows: 
Jq^ d(j) Jg p{9,(j),z',9\(p',X)sin{9)d9 = 1. This function gives the probabihty for a photon of wavelength 
A, incident on the scattering particle with angles {9', cf)') to be scattered in the direction {9, 4>) and satisfies 
p{9, (j>, 9' ,(/>', X) = p(cosG,A) where Q is the scattering angle, related to the polar and azimuthal angles by 
cos 9 — cos 9' cos 9 + sm9' sin0 cos(0 — (/)'). 

The solution of Eq. is splitted into two terms 

I{z, 9, cj), X) = I^'^iz, X) + fi, 0, A). (2) 

I'^''^{z, A), the direct intensity, is the solution of Eq. when there is no multiple scattering (no integral term): 

Mo^^ = -/3t(^)/''^'^(^) (3) 
where fiQ is the cosine of the polar angle for the incident radiation. The diffuse intensity is the solution of 

dI'''f(z,^^,<j>,X) ^ , f3sca{z,X) f'\^, 

— = -Pt{z,X)I A) + — J d(t> (4) 

dM'p(^,/i>,^',M',0',A)/'''^(z,/i',0') 

/3,ea(^,A)/*''(z,A) 
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-p{z,fi,(j>, z',-|^o|>o,A)', A). 



with fi = cos 6*. In this equation the variable /x takes values in the range — 1 < < 1. Negative ^ corresponds 
to radiation going downwards, whereas positive jj, describes radiation going upwards. In order to solve Eq. 
Q the diffuse intensity is expanded in a 2n Fourier cosine series (from now on we drop the superscript dif): 
I{z, ^,4>,X) = J2m=o -^'"(-z, m) cos m((/)o — if)). The phase function is expanded in a basis of 2n Legendre poly- 
nomials p(z, </), 2:', /z', (/)', A) = p(z,cos6,A) = ^^"q ^(2/ + l)5iPi(cos6). With these transformations and the 
theorem of addition of spherical harmonics Eq. Q becomes a set of integro-differential equations depending 
only on the z and fi coordinates: 

dl™-(z 11) 

^dT^ = -/3T(^,A)/"(z,/i) + J™(z,m), (5) 



where 

2n-l 

2 



j™(,,^) = ^££f^ ^ (2/ + l)5-p-(^) (6) 
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being P["ip) the associated Legendre polynomial, = gijfp^, and gi = \ Jj^^ p(cos 6)P;(cos Q)d(cos 6). 

In the discrete ordinate approximation the angular integral term in Eq. © is represented by a summation 
over n Gaussian quadrature points /i^ (fixed angles) also called " streams" . The intensity given by Eq. must 
satisfy boundary conditions at the top {z = z^ax) and bottom [z = 0) of the medium. Therefore we end up 
with a system of n coupled ordinary differential equations of the type ((HJ, one for each stream /Xg, and subject 
to a two-point boundary condition. Solving these equations, we will get a discrete approximation to the angular 
distribution of /'"(z, fj) from the top of the atmosphere to the surface level. 

In this work, for simplicity, we will consider the two-stream approximation to this problem and obtain the 
two streams (downwards) and /™(z,/i2) (upwards) fulfilling the boundary conditions that impose, 

first, no diffuse radiation incident at the top, and second, no radiation reflected back at the surface. This is 
expressed as follows, 

/"(^ = W,Mi) = (7) 
/"(z = 0,/i2) = 0. (8) 

In some cases it is preferred to use the optical depth of the medium r instead of the vertical geometric distance 
z. This is a non-dimensional variable which is defined as t(z) = PT{z)dz. It describes the attenuation within 
the medium of an incident beam of radiation when emission and multiple scattering are ignored, i.e. it is a 
measurement of the direct component attenuation. 

III. THE BENCHMARK PROBLEMS 

In this section we define three benchmark problems with an exact solution where the optical properties 
(3t{z, A) and (3sca{z, A) in Eq. ^ have different altitude (z) dependences. For the sake of simplicity we consider 
that particles scatter radiation uniformly, i.e. the scattering phase function is p(cos0) = 1 and that z^ax = 1- 

For the integral part in Eq. and following the two-stream approximation we consider a weighted summation 
over the streams fii ~ — 1/-\/(3) (down) and 112 = l/-\/(3) (up) which are the zeroes of the second order Legendre 
polynomial, i.e. P2{P'i) — 0, i = 1,2. Then, Eq. Q can be transformed in two coupled ordinary differential 
equations for the down (/i) and up {I2) radiation intensity: 

= + _ ^i, + __ j , (9) 

d7 ^ + — 2 — r ' ) ■ ^ ^ 

After a change of variables z' = I — z (from now on dropping the prime in z) and taking into account that 
Hi = —^2 the above equations may be rewritten as, 

Mo-^ = [A(z)-S(z)]Ai2/*'^ (11) 
az 

^ = A{z)h+B{z)l2+l'-{z)^ (12) 
az Ztt 

g = -A{z)l2 - Biz)h l'"-{z)^ (13) 



where 
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^^^^ ^ ~PT{Z)+Q.2bpsca{z) ^^^^ 

Bi^) = ^l^MfEfiM (15) 

r A(z)-B{z) , 

ld^-{z) = Ine-^ -0 ^^''^ (16) 

and /q is fixed by the condition at the top /''"'(z = 0). In the foUowing examples we will solve this problem with 
the two-point boundary conditions given in Eqs. Q"®, which are now explicitely written as /i(z = 0) = (at 
the top) and l2{z = 1) = (at the bottom). The incident intensity on the top is characterized by Iq = 100 and 
^0 — —0.788. We will consider two different height dependences of the optical properties. 



A. Linear dependence 



We propose as study-case a linear dependence on the extinction coefficients. In particular we have considered 
the functions A{z) and B{z), 

Aiz) ^ iii^ (17) 

B{z) = ^ ^ (18) 

With this choice the general solutions of Eqs. H12(l - (|13|l for the diffuse components are, 



16 7r(-c2 /i22 -h4a2 /zo 

for the downwards radiation, and 



2a — c 2a + c 16 7r(— ^2 +4 a^/^o^) 

for the upwards component. Here Ci and C2 are arbitrary constants fixed with the boundary conditions in Eqs. 
(0)-®. The direct component of the intensity is now 

/''^'■(z) = /oeV J . (21) 



For the specific choice c = —10.5, a = 4 the constants Ci and C2 that satisfy the boundary conditions are 
Ci — —33.0927 and C2 — — 6.71346e — 05. For this particular set of parameters the total optical path is r « 3. 
The optical functions (}t{z) and Psca{z) are shown in Fig. I21-I and in Fig. I21-II the attenuation of the direct 
intensity as a function of height. The two-stream components of the diffused intensity are shown in Fig. (SJ-I. 



B. Exponential dependence 

In atmospheric layers the gas density typically depends on the altitude as a growing exponential from z = 
(the top) to z = 1 (the bottom). For this reason the choice of exponential dependence in the extinction 
coefficients is very convenient as it allows us to validate our method with solutions similar to those appearing 
in atmospheres, where our method has been applied In particular we have chosen, 
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M^) - ^ ' (22) 

(az) / 2 52 _ 2N 

Biz) = ^ (23) 



The general solution for the downwards intensity is, 



Tt\ (e("-'&)^ , (-e(-)6)^ j-C^+a^b^) (/Xq^ -/Z0/Z2) ^ (- 

/i(z) = e^ ^Cg+e^ " T^— ^ — 2 2 , 2^2 2r-'o^^ 



(24) 



and for the upwards intensity is, 



c + ao — c + a6 16 tt (a^ 0^ /^o^ — ^2 ) 

where Ci and Ci are arbitrary constants fixed with the boundary conditions. The direct component of the 
intensity is, 

/*''-(z) = /oe("^^^). (26) 

For the specific choice c = —0.8, a = 3, 6 = 0.2, the constants C\ and C2 that satisfy the boundary conditions 
are C\ = —71.6787 and C2 — — 8.51812e — 05. The total optical path, for this particular set of parameters, is 
T w 3. The optical functions Pt{z) and Pt{z) are shown in Fig. and the direct intensity in Fig. The 
two-stream components of the diffused intensity are shown in Fig. 



IV. THE NUMERICAL METHOD 



In contrast to the examples discussed above, where the optical properties of the media are known functions, 
in the general case the medium optical properties are only known, based on measurements or models, at a 
discrete set of points. In order to solve the radiative transfer equation, it is generally assumed that these 
optical properties are piecewise constant in a layer around the point, i.e. that the plane parallel layers are 
homogeneous piecewise. As we shall see later this should not be assumed in the case of atmospheric layers, 
where the optical properties show a strong (exponential) altitude dependence. We therefore propose to solve the 
problem numerically, using an adaptative step procedure and interpolating the optical values where necessary. 

Equation (j^J is an initial value problem and it is easily integrated knowing the incident intensity on the 
top Iq. Equations (|12|I - H13|I and their boimdary conditions constitute a two-point boundary problem which 
can be solved for instance with a shooting method as those explained in |^. For problems with linear 
boundary conditions, shooting methods are able to get the exact solution in few steps. We choose values for all 
the variables at one boundary, which must be consistent with any boundary condition there, but are otherwise 
arranged to depend on arbitrary free parameters whose initial values are guessed. We then integrate the ordinary 
differential equations with initial value methods, arriving at the other boundary, and adjust the free parameters 
at the starting point that zeros the discrepancies at the other boundary. Thanks to this procedure the problem 
is reduced to the solution of an initial value problem where no assumption of piecewise homogeneity is required. 

The numerical integration of the equation for the direct component (j3Jl and of the set of n ordinary differential 
equations (O for the diffused streams is done with a fifth-order Runge-Kutta (RG) scheme with variable step 0] • 
Given a equation ^ — f(z, li) where Ii{z) is known, the RG scheme uses a weighted average of approximated 
values of f{z,Ii) at several points within the interval (z,z -I- dz) to obtain /^(z -I- dz). In contrast to other 
methods where piecewise homogeneity is required, this method takes into account the variation along the layer 
of integration of the diffused intensity streams, the direct component and the medium characteristics. For 




FIG. 3: Exact two-stream components Ii{z) (down) and hiz) (up) of the difTused intensity from the top z — to the 
bottom z = 1 of the medium when the scattering and absorbing functions are linear (I) and exponential (II). In both 
cases the total optical path is rjv « 3 and the incident direct radiation lo — 100. 



the adaptative step control we use a "step doubling" technique: the local error is estimated by comparing a 
solution obtained with a fifth-order scheme and the one obtained with a fourth-order method. The integrating 
step is halved if this error is above a desired tolerance. Next, given the tabulated function PT,i = Prizi) and 
Psca.i — Psca{zi) with 1 = I...N, wB interpolate to obtain the new values at the locations required by the RG 
and the adaptative step doubling integration scheme. We have chosen a cubic interpolating spline which allows 
the interpolating formula to be smooth in the first derivative, and continuous in the second derivative, both 
within the interval and at its boundaries |5|. 

As an example, next we solve for the direct and diffused intensity in the two benchmark problems explained 
in section ITTll For problems with optical depths of the order of r w 1, standard calculations based on the 
homogeneous layer approximation divide the medium into iV = 10 homogeneous layers j2]. We will analyze 
problems of r « 3 and keep a similar ratio for our comparisons, we will then increase the number of layers. The 
tabulated optical properties l3T{zi), (3sca{zi) at N equidistant Zi points are given as input. 

A summary of the maximal relative error E{I) — \I'2if"' — Idif'^^\l ^dif'^^ ^'^'^ t^*^ benchmark problems as a 
function of the method and number of initial layers is given in table We compare the relative error in the 
solution where piecewise homogeneity is assumed with the one obtained with our RG, step doubling, interpo- 
lating method. Figure 01 shows the relative error as a function of height z for the step doubling interpolating 
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technique with iV = 30 layers, for the hnear (I) and exponential case (II). For a ratio of layers N/t w 10 
the piecewise homogeneous approximation leads to errors of the order of 10% which can be reduced to 1% by 
increasing this ratio to N/t w 80. The step doubling interpolating RG technique permits the adaption of the 
numerical integration to the local characteristic length-scales of the medium and, for the same input, increases 
the accuracy of the solution between one and two orders of magnitude allowing for fast and accurate radiative 
transfer calculations. 



TABLE I: Solution of the linear and exponential problems with t « 3. Maximal local relative error in the diffused 
intensity as a function of the number of layers. 



linear case 


N=30 


N=60 


N=240 


piecewise homogeneous layers 


11% 


5.5% 


1.4% 


adaptative step interpolating method 


0.07% 


0.07% 


0.07% 


exponential case 


N=30 


N=60 


N=240 


piecewise homogeneous layers 


15% 


8% 


2% 


adaptative step interpolating method 


0.63% 


0.16% 


0.09% 



V. CONCLUSIONS AND DISCUSSION 



We have considered the discrete-ordinate approach to the radiation transfer equation. In view of our cal- 
culations the error in the plane parallel approach, assuming N piecewise homogeneous layers, with N « lOr 
can go up to 10% of the diffused intensity and thus the inhomogeneities in the atmosphere cannot be ignored. 
We have validated a general purpose numerical method that, based on an adaptative step integration and an 
interpolation of the local optical properties, can significantly improve (up to two orders of magnitude) the ac- 
curacy of the solution. This is furthermore of interest for practical applications, such as atmospheric radiation 
transfer, where the scattering and absorbing properties of the media are only known, based on experiments 
or theoretical models, at certain discrete points. Furthermore, this numerical method can be straightforward 
extended to multiple streams, non trivial boundary conditions and non uniform scattering function allowing for 
fast and accurate solutions of the radiative transfer equation. 
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